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SUMMARY 

This report documents the changes that were made to the two-equation k-e turbulence model 
in the NPARC (National-PARC) code. The previous model, based on the low Reynolds number 
model of Speziale, was replaced with the low Reynolds number k-E model of Chien. The most 
significant difference was in the turbulent Prandtl numbers appearing in the diffusion terms of the k 
and E transport equations. A new inflow boundary condition and stability enhancements were also 
implemented into the turbulence model within NPARC. The report provides the rationale for making 
the change to the Chien model, code modifications required, and comparisons of the performances of 
the new model with the previous k-£ model and algebraic models used most often in PARC/NPARC. 
The comparisons show that the Chien k-£ model installed here improves the capability of NPARC to 
calculate turbulent flows. 
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SYMBOLS 

constant for turbulent kinetic energy dissipation expression (equal to 0.164) 

skin friction along flat plate 

k-£ turbulence model terms 

k-£ turbulence model constant (equal to 0.09) 

k-£ turbulence model terms 

distance from centerline to top or bottom wall of ejector nozzle 
turbulence intensity 
turbulent kinetic energy 
turbulent mixing length 

Reynolds number based on turbulent quantities 
Reynolds number based on momentum thickness 
time 
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U, V, w 

velocities 

t 1 t 

U , V , w 

turbulent velocity components 

u ref(k-£) 

reference velocity for turbulent kinetic energy limiter 

Utot 

local inflow velocity magnitude 

x > y 

Cartesian coordinates 

y+ 

distance from wall normalized by shear length scale 

e 

rate of turbulent kinetic energy dissipation 

P 

dynamic viscosity 

Pt 

turbulent viscosity 

Mr- max 

maximum turbulent viscosity limiter 

n 

production term in k-e model 

p 

density 


turbulent Prandtl numbers 

X 

turbulent time scale 

Subscripts: 

ij 

computational coordinates 


INTRODUCTION 

The PARC Navier-Stokes code (refs. 1 and 2) has been used extensively by government and 
industry to analyze propulsion flows. Despite advances in flow-solving capabilities, the ability of 
codes such as PARC to calculate complex flows is strongly dependent on the turbulence model 
employed. Until recently, only algebraic models have been available in PARC for turbulent flow 
simulations. The standard algebraic turbulence model in PARC is based on the work of P.D. Thomas 
(ref. 3). This model calculates turbulent viscosity near surfaces (wall-bounded part of the model), 
and in regions where two or more flows are mixing (free shear layer part of the model) but was 
optimized for the latter. The Baldwin-Lomax algebraic turbulence model (ref. 4), which only 
calculates turbulent viscosity in wall-bounded regions, is available also. Algebraic turbulence models 
such as these often model complex flow cases inadequately because the single mixing length 
distributions used to calculate turbulent viscosity are normally tuned to a particular case and often are 
not applicable to all flows. 

Two-equation models avoid this single mixing length limitation by solving two additional 
transport equations to calculate turbulent viscosity. These models have been installed previously in 
PARC to improve the code’s capability to calculate turbulent flows. The Chien low Reynolds number 
k-e model (ref. 5) with modifications for compressibility added by Nichols (ref. 6) has been available 
in the 2D/axisymmetric code (PARC2D) but was not successfully installed in PARC3D. Another k-e 
model based on the work of Speziale (ref. 7) was installed in PARC2D and PARC3D but has not 
provided desirable accuracy or stability. 

To provide the U.S. aerospace community with a reliable Navier-Stokes propulsion flow 
simulator, the NPARC (National-PARC) Alliance was recently instituted and the PARC code has been 
renamed to NPARC. As part of this effort, modifications to the two-equation turbulence model in 
NPARC were made to improve the model’s accuracy and to enhance its stability. This report 
documents these modifications. The following sections provide the rationale for making the required 
code modifications, recommendations for using the model, and comparisons of the new model with 
the turbulence models most often used in PARC/NPARC for two benchmark cases. 


BACKGROUND OF SPEZIALE AND CHIEN k-e TURBULENCE MODELS 

The two-equation turbulence model in the NPARC code (Version 1 .0) was installed by 
Nichols (ref. 8) in 1991. It is a k-e model (two equations are solved: one for turbulent kinetic 
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energy k and the other for the rate of turbulent kinetic energy dissipation e) based on the work of 
Speziale (ref. 7). In reference 7, Speziale develops a k-x two-equation turbulence model (the 
quantity solved in the second equation, x, is a turbulent time scale and is equal to k/e) to avoid solving 
for e, which he states does not have a natural boundary condition near solid surfaces. Speziale, 
however, shows that his x-transport equation is equivalent to the e-transport equation if an assumption 
holds (to be discussed shortly). This provides the following k-e model, which has been available for 
use in NPARC: 

pk 2 /e (1) 


D(pk) 3 ( p t 1 3k 0 

VK - ~ = — |i + — + n - pe 

Dt dXi 1 cjaxj 


= — I + c £l f,n^- -c £ 2 f 2 p^ 

k k 


Dt 9Xj ^ o e J 3x ; 


3u: 5uj 9u. 

n = — 1 + — i 

dxj 3x ; 3xj 


f = (1 + 3.45 / ^/Re t ) tanh (y + / 70) 


f, = 1.0 



(7) 

( 8 ) 


with C„ = 0.09, C el = 1 .44, C e2 = 1.83 ( 1 - .22 e(-Re, /6)2), a k = l .36, and 0 e = 1 .36^ 

The constants 0 k and o e are turbulent Prandtl numbers appearing in the diffusion terms ot 
the k and e transport equations (2) and (3). Speziale’ s assumption allowing his second transport 
equation for x to be transformed to an e transport equation is that o k = o e = 0 t Launder and 
Spalding (ref. 9) and Launder et al. (ref. 10) indicate that the values for the turbulent Prandtl 
numbers 0 k = 1.0 and o e = 1.3 are appropriate for flows involving mixing layers. Reference 9 
mentions that values other than these may have been used where strictly wall-bounded flows were to 
be calculated (as was Speziale’ s intent) but that 0 k = 1.0 and o e = 1.3 are more applicable to a wide 
class of flows involving wall-bounded and/or mixing regions. From equation (1) and the diffusion 
terms of equations (2) and (3), note that for a given flow field, increasing o k substantially (from 1.0 
to 1.36) while holding 0 e to a smaller change (from 1.3 to 1.36) results in lower turbulent viscosity p t , 
and possibly significantly less calculated mixing for a case involving mixing of flows. Because the 
low Reynolds number k-e model of Chien (ref. 5) employs 0 k = 1.0 and 0 e = 1.3 and has been used 
successfully in many Navier-Stokes codes, the current work incorporates the Chien k-e model in 
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place of the existing k-e model in NPARC. For reference, the k-e model described in equations (1) 
through (8) will be referred to as the NPARC 1 .0 k-e model in the rest of this report. 

The Chien low Reynolds number k-e model was installed successfully in the PARC2D code in 
1990 (refs. 6 and 11). It is presented next to demonstrate the changes required to convert the 
NPARC 1.0 k-e model to the Chien model in NPARC: 
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with C^ = 0.09, C e , = 1.35, C e2 = 1.8, c k = 1.0, and a £ = 1.3. 

By comparing equations (1) through (8) to equations (9) through (16), note that the Chien 
and NPARC 1.0 k-e models differ in C el , C e2 , o k , o e , f^, and f 2 . The Chien model also has an extra 
near-wall term for each of the two turbulent transport equations (10) and (11) that does not appear in 
the NPARC 1.0 k-e model. These terms allow for all of the turbulent quantities to be set to zero at a 


solid surface by using the Chien model. These changes were made to convert the k-e model in 
NPARC from one based on the Speziale model to one based on the Chien model. 


INFLOW BOUNDARY CONDITIONS 
Free Inflows 

The standard boundary condition in NPARC for subsonic inflows is the free boundary 
(type 0). The previous implementations of the Chien and NPARC 1.0 k-e turbulence models 
extrapolated the quantities k and e from the interior when this free boundary was specified at an 
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inflow. This has been shown to result in artificially low levels of turbulent kinetic energy and 
calculated turbulent viscosity near inflows (ref. 12). 

To allow for more accurate simulation of turbulence at an inflow when using the free 
boundary, a new inflow boundary condition for the k and £ equations was added in the current work. 
Turbulent kinetic energy is defined as 


If isotropic turbulence is assumed, 

I 2 

so that the turbulent kinetic energy may 

k = ^I 2 • ufo. (19) 

The implementation of this boundary condition in NPARC requires the intensity I to be specified 
which allows k to vary as the local velocity magnitude at the free boundary inflow. With k obtained 
from equation (19), one possibility for specifying E is to set the turbulent viscosity equal to the 
laminar viscosity or a multiple of the laminar viscosity, as suggested in reference 13. Equation (9) is 
then used to determine the corresponding value of £ at the boundary. A second method for 
obtaining £ is to use the following expression from reference 14, which requires a turbulent mixing 
length l to be specified: 

£ = C D k 3/2 /* (20) 

The input of the turbulent intensity and either the turbulent viscosity or turbulent length scale 
required when using this free boundary condition is described in the code usage section. 

Fixed Inflows 

If a fixed boundary (type -10), which is recommended for supersonic inflows, is specified at 
an inflow, all main flow and turbulent quantities remain unchanged from the values on the restart file. 
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( 18 ) 


be defined as a function of the turbulence intensity: 


STABILITY ENHANCEMENTS 

Experience with the NPARC code has shown that the previous k-£ model had stability 
difficulties, especially in the NPARC3D code. To enhance the stability of the current k-£ model, the 
following modifications were made to the NPARC code: 

1. With the k-£ model previously installed in NPARC, only the turbulent viscosity was 
relaxed, and for 50 iterations following initialization of the k-£ model from an algebraic turbulence 
model. ’ In the current Chien model, the updated turbulent quantities k and £ were relaxed in addition 
to the turbulent viscosity, and for 500 iterations after initialization of the k-£ model. This was the 
method used successfully in PARC2D (ref. 11). 

2. Examination of unstable flow cases indicated that the instability was caused by 
unrealistically large values of k which often were calculated during the initial iterations after switching 
from an algebraic model. Examining equation (9), very large values of k cause p t to become very 
large. To prevent such large values of k and resulting code was added to limit any value of k from 
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exceeding 10 percent of a reference kinetic energy obtained from a velocity 

( u ref(k-e)) : 


k = minimum ^k - calculated; 0.10 • ~u^ f(k _ £) J 


characteristic of the flow 

( 21 ) 


In addition, the turbulent viscosity is prevented from exceeding a maximum turbulent viscosity 
1-h-max- Input of this characteristic velocity and maximum turbulent viscosity is discussed in the next 
section. 


CODE USAGE 

The NPARC code with the Chien k-e model implementation is run essentially as described in 
reference 8. The k-e model is initialized from an algebraic turbulence model solution at the iteration 
level set by the input NTURB. As previously mentioned, all turbulent quantities are relaxed for the 
first 500 iterations after NTURB to assist stability. The limiters in the code are set with the inputs 
UREFKE and TMUMAX. UREFKE is the characteristic velocity nondimensionalized by AREF (the 
reference speed of sound for NPARC), and TMUMAX is the maximum turbulent viscosity 
nondimensionalized by the reference viscosity. Default values for these inputs are UREFKE =1.0 
and TMUMAX = 3000.0. These should be reasonable upper bounds for many flows of interest; 
however, if the flow to be calculated is at very high or very low Mach number, the UREFKE variable 
can be set to a velocity considered characteristic of the flow to be calculated. In addition, if the flow 
contains a high-speed free shear layer, TMUMAX can be increased appropriately. 

The new inflow boundary condition (used with the main flow solver’s free boundary) allows 
up to 3 different inflows to be set in namelist TURBIN. Values for intensity are specified through the 
inputs TUIN1, TUIN2, TUIN3, and the corresponding turbulent viscosities can be set through the 
inputs TMUIN 1 , TMUIN2, and TMUIN3. For example, a flow has two inflow boundaries. The first 
has a measured turbulence intensity of 2 percent and the turbulent viscosity is to be set equal to the 
laminar viscosity; the second has a turbulence intensity of 5 percent and the turbulent viscosity is to 
be set to 10 times the laminar viscosity. The inputs would be set as follows: 

&TURBIN 
TUIN1 = .02, 

TUIN2 = .05, 

TMUIN 1 = 1., 

TMUIN2 = 10., 

&END 

Also, the EDGE (JEDGE, KEDGE, or LEDGE) parameters must be set to -1 for the first boundary 
and -2 for the second boundary. A value of -3 would have been specified if there were a third inflow 
boundary. Negative integers are used to avoid confusion with the EDGE used for no-slip surfaces. If 
anything other than -1, -2, or -3 is specified for an EDGE input, the default extrapolation boundary 
condition is used for k and £ at that inflow boundary. A turbulent length scale can be specified 
instead of a turbulent viscosity by specifying the length (nondimensionalized by the NPARC 
reference length) as a negative value with the TMUIN1, TMUIN2, or TMUIN3 input. For example, at 
an inflow boundary, the turbulence intensity is 2 percent and the desired turbulent mixing length is 
0.05 ft (with the reference length equal to 1 ft). The inputs would be 

&TURBIN 
TUIN1 = 0.02 
TMUIN 1 = -0.05 
&END 
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If a fixed boundary (type -10) is specified, the turbulent quantities will remain the same as those 
values on the restart file at the boundary, just as for the main flow quantities. Outflow boundary 
conditions can be specified with the free boundary (type 0) and nothing set for the EDGE 
parameters. The free boundary condition will extrapolate the turbulent quantities, which is 
appropriate at outflows. Slip surfaces or planes of symmetry can be modeled with the slip wall 
boundary (type 50) and axes of symmetry are modeled with the type 51 boundary in the NPARC2D 
code. These boundary conditions also extrapolate the turbulent quantities from the flow interior. 
The no-slip boundary conditions (types 60 and 61) set all of the turbulent quantities to zero at the 
boundary. The pole boundaries (types 81 to 83) and block interface boundaries (types 70 to 75) 
operate the same for the turbulent flow quantities as for the main flow quantities. 


RECOMMENDATIONS 

Four recommendations for using the Chien k-£ model are (1) All INVISC values, which 
control the directions for which viscous terms are calculated, should be set to 1. (2) When the k-£ 
model is initialized, DTCAP should not be increased from the value that was locked on while using 
the algebraic model. The time step should be examined after turning on the Chien model, especially 
during the iterations immediately after the iteration level corresponding to (NTURB + 500) as it may 
change from what was acceptable when using the algebraic model. DTCAP will most likely need to 
be reduced slightly. (3) The grid should be packed to solid walls such that the first grid point off of 
the wall is within y+ = 5 and preferably positioned closer to y+ = 1. Avva, Smith, and Singhal 
(ref. 15) have shown the sensitivity of low Reynolds number k-£ models to the grid packing near 
solid surfaces and the importance of locating the first point away from the wall well within the 
laminar sublayer. (4) As a calculation approaches convergence, the flow field should be examined to 
be sure that UREFKE and TMUMAX are not limiting the turbulent kinetic energy and turbulent 
viscosity. These inputs are intended only to assist stability of the Chien k-£ model in the iterations 
after initialization from an algebraic turbulence model. 

It is necessary to run all of the blocks in turbulent mode when using the k-£ model because it 
is not possible to interpolate the turbulent quantities from an inviscid or laminar block. 


INVESTIGATION OF NPARC TURBULENCE MODELS 

Two benchmark test cases are investigated in this section to compare the performance of this 
Chien k-£ turbulence model, the previous NPARC 1.0 k-£ model, and the algebraic turbulence 
models used most frequently in PARC and NPARC. The first case is flow over a flat plate to 
investigate the models’ performance for a wall-bounded flow, and the second is flow through an 
ejector nozzle to investigate a more complex propulsion flow case with the mixing of two flows as the 
dominant flow feature. Both flows are two-dimensional so that they could be properly modeled with 
NPARC2D or NPARC3D. The Baldwin-Lomax model was not applied to the ejector nozzle case 
because it is a wall-bounded turbulence model and this flow is dominated by turbulent mixing. 
Calculations were obtained with both versions of the code for each test case and turbulence model. 
However, because the results obtained by using NPARC2D or NPARC3D for each case were very 
similar, only the NPARC3D results are presented in the following discussion. 

Flat Plate 

A schematic of the flat plate case investigated with NPARC is shown in figure 1. The inflow 
total pressure and outflow static pressure, both free boundaries, provided a Mach 0.2 flow over the 
plate. The grid had 111 points in the horizontal direction and 81 points in the vertical direction and 
was packed tightly to the wall to resolve the boundary layer. Five points were used in the z-direction 
to accommodate the NPARC3D code. A fixed inflow profile was not used. Instead, the flow reached 
the leading edge of the plate at grid point 15 in the horizontal direction where the boundary layer 
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began and then continued downstream over the no-slip surface. Two inflow boundary conditions 
were specified for use with the new Chien model: (1) the default extrapolation, and (2) turbulence 
intensity set to 5 percent and the turbulent viscosity set to 100 times the laminar viscosity. This 
second inflow was specified to determine if setting a relatively high turbulence intensity (5 percent) 
would have an effect on the flow calculation, and not because these specific turbulence quantities 
were measured in an experiment. The calculations showed that the flow over the flat plate was 
independent of the upstream condition. 

NPARC3D calculations are compared with experimental data of Wieghardt (ref. 16) in 
figure 2. The skin friction calculated by using the Thomas model is very low because of the very 
small turbulent viscosity values that the Thomas model calculates for purely wall-bounded flows 
(such as this flat plate test case). The Baldwin-Lomax model produces very good agreement with 
experimental data as expected because it was tuned to attached wall boundary layer flows. The 
NPARC 1 .0 k-e model provides reasonable agreement with the data and the Chien k-e model 
(installed in this work) provides somewhat closer agreement. The Chien solution shown in figure 2 
was obtained by using the extrapolation inflow boundary condition for the turbulent quantities, but 
both Chien calculations provided the same flow solution as mentioned earlier. 

Ejector Nozzle 

A schematic of the ejector nozzle test case is shown in figure 3. This flow is dominated by 
mixing of the primary flow with entrained secondary flow and is similar in fundamental operation to 
nozzles that have been considered for vertical and short takeoff or landing (VSTOL) and high-speed 
transport application. The grid had 131 points in the horizontal direction by 121 points in the 
vertical direction. Five points were used in the z-direction with NPARC3D, as for the flat plate case. 
Free boundaries were specified at the inflow to both the primary nozzle and ejector inlets and at the 
outflow. Initial calculations with the NPARC 1 .0 k-e and new Chien k-e models used the default 
extrapolation boundary because no turbulence measurements were available from the experiment. 
Atmospheric pressure was specified at the ejector inflow while the total pressure specified for the 
primary nozzle was set to provide a nozzle pressure ratio (nozzle total pressure divided by 
atmospheric static pressure) of 2.44. The outflow pressure was the static pressure measured in the 
experiment (ref. 17) at approximately 26.9 cm downstream of the primary nozzle exit. Because the 
geometry of the nozzle is symmetric about a plane passing through the center of the ejector nozzle 
assembly, only one half of the ejector nozzle was modeled with NPARC. 

Velocity profiles obtained from the NPARC calculations are compared to experimental data 
in figure 4. The axial positions were measured relative to the primary nozzle exit plane and the 
vertical positions were nondimensionalized by the local half-duct height H. The velocity profiles at 
7.6, 12.7, 17.8, and 26.7 cm downstream of the primary nozzle exit plane show that the Chien k-e 
model added in this work provided much better agreement with the experimental data than that 
provided by the Thomas algebraic model (which was optimized for free shear layers such as in this 
ejector nozzle case) and the NPARC 1.0 k-e model, despite the inflow boundary condition. Total 
temperature profiles in figure 5 indicate the same mixing behavior. The Thomas and NPARC 1.0 k-e 
models produce significantly less turbulent viscosity than the Chien model (fig. 6) which is directly 
responsible for the large differences in predicted mixing. The large differences in the two k-e 
solutions for this mixing-dominated problem is due to the specification of constants and o e . As 
mentioned earlier, the values for the turbulent Prandtl numbers used by the Chien model, = 1 .0 
and c e = 1.3, have been shown to provide more realistic flow predictions than the values = 1.36 
and g £ = 1.36 used by the NPARC 1.0 k-e model. The enhancements added in this work gave the 
Chien k-e model more stability than the NPARC 1 .0 k-e model in the iterations immediately 
following initialization from the algebraic model. 

Additional calculations were obtained with the Chien k-e model while varying the inflow 
boundary condition to determine the effects of specifying the turbulent quantities on the mixing 
downstream. In addition to the extrapolation case (the default inflow) previously discussed, two other 
inflows were examined: (1) turbulence intensity = 5 percent and turbulent viscosity = 500 times the 
laminar viscosity for the primary flow, turbulence intensity = 2 percent and turbulent viscosity = 100 
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times the laminar viscosity for the secondary flow; and (2) turbulence intensity = 5 percent and a 
turbulent length scale set to 5 percent of the nozzle exit height for the primary flow, turbulence 
intensity = 2 percent and turbulent length scale = 5 percent of the secondary inlet passage height at 
the primary nozzle exit plane. These two additional inflows were specified arbitrarily, since no 
turbulence measurements were available from reference 17. 

Figures 7 and 8 compare the solutions obtained with the Chien k-e model and the three 
(jjffergjit inflow boundary conditions. The calculation that specified the turbulence intensity and 
turbulent viscosity for the primary and secondary flows more closely matches the data than the 
calculation using extrapolation at the inflows or the calculation specifying the turbulence intensity 
and turbulent mixing length. The specified turbulence intensities and viscosities cannot be justified 
based on these results because these quantities were not taken from the experiment. However, these 
results do show that the new inflow boundary condition may enable more accurate calculations for 
such flow cases when turbulence measurements are available to set as inflow boundary conditions. 


CONCLUSIONS 

The two-equation turbulence model in the NPARC code was modified so that the model is 
based on the low Reynolds number k-e model of Chien and no longer is based on the Speziale 
model. The most significant change was made to the turbulent Prandtl numbers appearing in the 
transport equations for k and e. Stability enhancements and a new inflow boundary condition for the 
turbulent quantities were also added to the k-£ model. Comparisons of the NPARC solutions 
obtained using the previous and new k-£ models with experimental data indicate that the Chien k-£ 
model installed in this work improves the capability of NPARC to calculate propulsion flows. In 
addition, the limiters added for code stability seem to improve the convergence characteristics of 
NPARC relative to that observed with the NPARC 1.0 k-£ model. The specification of different 
inflow boundary conditions resulted in different mixing solutions for the ejector nozzle case. The 
new inflow boundary condition will be examined further to establish guidelines for appropriate 
specification of the turbulent quantities at inflows. 
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Figure 1. Schematic of the flat plate test case. 
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Figure 2. Skin friction for the flat plate flow. 
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Figure 6. Turbulent viscosity contours for the 2D ejector nozzle. 
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(a) x = 7.6 cm. (b) x = 26.7 cm. 

Figure 8. Total temperature profiles for the 2D ejector nozzle flow 
using different inflow boundary conditions. 
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